Quantitative test of thermal field theory for Bose-Einstein condensates 
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We present numerical results from a full second order quantum field theory of Bose-Einstein 
condensates applied to the 1997 JILA experiment [D. S. Jin et al, Phys. Rev. Lett. 78, 764 
(1997)]. Good agreement is found for the energies and decay rates for both the lowest-energy m — 2 
and m = modes. The anomalous behaviour of the m — mode is due to experimental perturbation 
of the non-condensate. The theory includes the coupled dynamics of the condensate and thermal 
cloud, the anomalous pair average and all relevant finite size effects. 
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One of the most intriguing consequences of the experi- 
mental realization of Bose-Einstein condensation (BEC) 
was the prospect of quantitative tests of finite tempera- 
ture quantum field theory (QFT). The pioneering mea- 
surements of condensate excitations at JILAprovide the 
most stringent tests to date of such theories 0, y] . Accu- 
rate calculations have proved difficult, however, because 
of the need to include the dynamic coupling of condensed 
and uncondensed atoms simultaneously with effects due 
to strong interactions and finite size. In this paper we 
describe the first direct comparison of a full second or- 
der QFT calculation with the JILA measurements. Our 
results show that accurate tests of QFT are possible if it 
is properly adapted to the finite, driven systems under 
consideration. 

Measurements of excitations at low-temperature are 
in good agreement with predictions based on the Gross- 
Pitaevskii equation (GPE) and Bogoliubov quasiparticles 
UlEl El' However, the finite-temperature JILA results Q 
have proved much harder to explain. In this experiment 
the energies of the lowest-energy modes with axial angu- 
lar momentum quantum numbers m — 2 and m — were 
measured as a function of reduced temperature t — T/T®, 
where T is the absolute temperature and Xj? is the BEC 
critical temperature for an ideal gas. The m = 2 mode 
was observed to shift downwards with t, while the m = 
mode underwent a sharp increase in energy at t ~ 0.6 
towards the result expected in the non-interacting limit. 

The temperature dependence of the excitations has 
been studied theoretically using the Popov approxima- 
tion to the Hartree-Fock-Bogoliubov formalism, where 
the anomalous (pair) average of two condensate atoms is 
neglected. This gives good agreement with experiment 
for low temperatures but can not explain the results for 
t > 0.6 [3- Good agreement for all t for the m = 2 mode 
was obtained using an extension of this approach which 
includes the anomalous average [(|, and also within the 
dielectric formalism Q. However, both approaches were 
unable to explain the upward shift of the m = mode, 
and an analytical treatment of the problem also predicted 



downward shifts for both modes pj. The importance of 
the relative phase of condensate and non-condensate fluc- 
tuations was emphasized by Bijlsma, Al Khawaja and 
Stoof (BKS) [9|, who showed that the experimental re- 
sults for m = can be qualitatively explained by a shift 
from out-of-phase to in-phase oscillations at high tem- 
perature. Jackson and Zaremba (JZ) obtained good 
quantitative agreement for both modes using a GPE for 
the condensate coupled to a non-condensate modelled by 
a Boltzmann equation. However, this approach neglects 
the phonon character of low-energy excitations as well 
as the anomalous average and Beliaev processes. The 
anomalous average can be significant, especially near a 
Feshbach resonance 0, 0] and Beliaev processes have 
been directly observed in a number of recent experiments 
|TSIT3.IT^ . It is therefore important to explain the JILA 
results using a theory which includes these effects. 

In this paper we present numerical results for the exci- 
tations of a dilute gas BEC at finite temperature for the 
conditions of the 1997 JILA experiment [2(. We find good 
agreement with the experimental results for both the 
m = 2 and m = modes, and in particular we are able 
to explain straightforwardly the anomalous behaviour of 
the m = mode. The results are based on a theoretical 
treatment recently developed by one of us (SM), as an 
extension of an earlier second-order perturbative calcu- 
lation [HI Hfl | . The formalism adapts the linear response 
treatment of Giorgini and closely models the experi- 
mental procedure where excitations are created by small 
modulations of the trap frequencies. The result is a gap- 
less extension of the Bogoliubov theory which includes 
the dynamic coupling between the condensate and non- 
condensate, all relevant Beliaev and Landau processes 
and the anomalous average. It is also consistent with 
the generalized Kohn theorem. The theory is valid in 
the collisionless limit of well-defined quasiparticles. For 
homogeneous systems at finite temperature this requires 
(fc B T/n [/o)(noa 3 ) 1/2 < 1, where uq is the condensate 
density, a s is the s-wave scattering length, fee is Boltz- 
mann's constant and Uq = AiTh 2 a s /m where m is the 
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atomic mass 0, ^J. For the JILA experiment Q this 
parameter does not exceed 0.03 at the trap centre for the 
highest temperature we consider. 

The theory starts from the generalized GPE for the 
condensate wavefunction <E>(r, t) 

'" P \H sp + P(r, t) - A(i) + N (t)U m 2] 



ih- 



dt 



+2U n{r, t)& + U m(r, - /(r, t). (1) 

Here H sp = —h 2 \7 2 /2m + Vt rap (r) is the static single- 
particle Hamiltonian, P(r, t) is the time-dependent exter- 
nal perturbation and X(t) is a scalar. The non-condensate 
density n(r, t), anomalous average m(r, t) and /(r, t) 
are constructed from time-dependent quasiparticle wave- 
functions Ui(r,t) and Vi(r,t) by 



n(r, £) 


= ^ Mr,t)| 2 ^ + |«i(r, i)| 2 (A^ + 1), 

i 


(2) 


m(r, t) 


= 5^u i (r,tK(r,t)(22V i + l), 


(3) 


/(r,*) 


= £c?(t)JVjUj(r,i) +c i (t)(AT i + lK(r,i), 


(4) 


Ci(t) 


= (7 /dr|*| a [#*u i (r,*) + *« i (r ) t)] . 


(5) 



The quasiparticle wavefunctions evolve according to 
d f Ui \ ( L M 



ih- 



dt \v 



-M* -U 



L(r, t) = H sp + P(r, t) + N U [|$| 2 + Q|$| 2 Q 
M(r,f) = AWoQ* 2 Q*, 



(6) 

(7) 
(8) 



where the orthogonal projector Q(r, r',t) = S(r — r') — 
4>(r, t) Jdr' $*(r', t) ensures orthogonality of the conden- 
sate and non-condensate. 

The quasiparticle populations {iVj} are independent 
of time and given by the Bose-Einstein distribution 
Ni = l/(e e, / feBT — 1) where e, is the Bogoliubov en- 
ergy (see below). Most quantities in the theory depend 
on temperature via these populations. The condensate 
population No(t) is defined in terms of the fixed to- 
tal number of particles N by No(t) = N — Jdrn(r,t). 
The zero-temperature part of the anomalous average 
rh(r,t) is ultra-violet divergent, but it can be renormal- 
ized straightforwardly [8Hli|. 

The above equations are obtained using the number- 
conserving approach of Castin and Dum, modified for 
finite temperature calculations [IlGlt. The terms f(r,t) 
and Q are a feature of this approach and do not appear 
in symmetry-breaking theories. We find that they can 
give a significant contribution to the energy shifts. 

In the static case, Eq. Q has a time-independent so- 
lution $(r,t) = $(r) which satisfies 



H, 



Bp -\ + N U \${r)\ 2 +2U h(r) *(r) (9) 



where A is the condensate eigenvalue, roughly equal to 
the chemical potential. If we set n, fh and / to zero, 
we obtain the usual GPE, with wavefunction $o(r) and 
energy Ao We solve Eq. by linearizing the change in 
energy and shape relative to this solution. Writing $ — > 
$o( r ) in Eq. ©, we obtain static quasiparticle wavefunc- 
tions Uj(r, t) = Uj(r)e~" lt / fi ', u<(r,t) = Vi{r)e~ leit / h and 
the Bogoliubov energies {&;}. These solutions are used 
to construct n(r), m(r) and /(r) and also provide a con- 
venient basis for the subsequent calculation. 

The external perturbation P(r, t) leads to all quanti- 
ties developing a small time-dependent oscillation around 
their static values, $(r, t) = <f>(r) +<5$(r, t), n(r, t) = 
h(r) + Sh(r, t), etc. Substituting these expressions into 
Eq. (JJJ and linearizing, we obtain the equation of motion 
for the condensate fluctuation <5$(r, t). This equation 
can be solved by combining it with its complex conju- 
gate, Fourier transforming and expanding the fluctuation 
in the static quasiparticle basis 



<5$*(r, -uj] 



Ui{v) 
Vi{v) 



(10) 



The expansion coefficients fcj(u>) are directly related to 
the condensate density fluctuations Sno = S(Nq\^ | 2 ), 
which are measured experimentally. 

Dynamics of the non-condensate can occur via two dis- 
tinct mechanisms; either it is driven directly by the per- 
turbation or indirectly via the condensate. If we neglect 
the first possibility and assume that only the single mode 
'p' is excited, then b p {ui) is given by 



b p (u) = P p o(uj)Gp(ll> + i-i). 



(11) 



Here P p q(lo) is the matrix element for the generation of 
the excitation from the condensate and 17 is a small imag- 
inary part in the frequency (discussed below) . The resol- 
vent G p (lu) is defined in terms of a frequency-dependent 
self-energy by 



E p (lj) 



huj — e p — S p (w) ' 



(12) 
(13) 



+C/ m(r)$*(r)-/(r) = 0, 



S p (w) contains two types of energy shifts, static (S) 
and dynamic (-D), corresponding to the different roles 
of the thermal cloud. The static term AE p S ^ comes from 
interactions between a condensate fluctuation and the 
static non-condensate mean-fields. The dynamic term 
AE p (ui) describes the driving of non-condensate fluc- 
tuations by the condensate and their subsequent back 
action, which leads to damping and energy shifts of con- 
densate excitations. The inclusion of this contribution 
leads to a gapless excitation spectrum 0, 0] . 

However, the non-condensate can also be excited di- 
rectly by the external perturbation, and can then gen- 
erate condensate excitations. This process therefore 
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changes the effective excitation matrix element P p q and 
can be included by replacing Q v (uj + if) in Eq. Ijllfl with 
the modified resolvent 1Z P (lu + ij) , defined by 



pO 



P (w). (14) 



The important extra term here is AP^'(tJ) which de- 
scribes the generation of non-condensate fluctuations by 
the perturbation and their subsequent coupling to the 
condensate A p ^ 
the static condensate shape [$0(1") 



AP p Q (w) describes the effect of changes in 



$(r)] 



The detailed definition of AE p D ^ and APp 1 ^' is lengthy 
and is given elsewhere 0, 0] . We note here that they 
are both calculated as a sum over many Landau and 
Beliaev processes which are resonant whenever an en- 
ergy matching criterion is satisfied. The parameter 7 in 
Eq. Ijlll) is required to keep AE p D ^ and AP p D ^ finite at 
the resulting poles. Its inclusion can be formally justified 
from the finite experimental resolution and its value is of 
order the inverse of the experimental observation time. 
Our numerical results are essentially independent of this 
parameter for physically relevant values. 

If E p and AP p D ^ are roughly independent of frequency, 
the energy shift can simply be calculated from the poles 
of Q, i.e. the solutions to E p = frw p = Real [e p + T, p (u> p )], 
while the decay rate is given by T p = — Imag [T, p (lu p )] jh. 
This situation arises in homogeneous condensates where 
an excitation couples to a continuum of decay channels 
and the resolvents are Lorentzians. For a finite system, 
however, T, p (co) depends on frequency, and neither Q p (ui) 
nor lZ p (cu) are perfect Lorentzians. In this case, we ex- 
tract energies and decay rates by fitting b p (uj) to a com- 
plex Lorentzian plus a constant (7 is subtracted from 
the resulting decay rate) . This corresponds to the exper- 
imental procedure of fitting a decaying sinusoid to the 
condensate density fluctuations in the time domain. The 
frequency dependence of P p q{lo) is included as a (known) 
weight function in the fit to ensure that only the experi- 
mentally relevant range of frequencies is included. 

We present numerical results for the parameters of 
the JILA experiment Q- We consider a condensate of 
6000 Rb atoms in an anisotropic harmonic trap with 
radial and axial trap frequencies of uo r = 2tt x 129Hz, 
u> z = 2ir x 365Hz. The scattering length is a s = 110 
Bohr. The condensate population is fixed for all the 
temperatures considered, which is consistent with the 
experimental results for t < 0.9. Zero-temperature ef- 
fects have been included using the appropriate ultra- 
violet renormalization. The external perturbation has 
the form P(r, t) cx r 2 cos(m p — a;dt)©(t)0(T<z — i) where 
r and <f> are the radial and azimuthal angle coordinates, 
~ e p /fi is the central drive frequency, O(r) is the unit 
step function and Td = 14ms is the drive time. The pa- 
rameter 7 is taken to be 7 = 0.036So; r 0). 
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FIG. 1: Ab initio theoretical excitation energies E (open sym- 
bols) compared with experiment (filled circles) for (a) m v — 



and (b) m v 



2. Diamonds neglect direct thermal driv- 



ing (G P ), open circles include it (1Z P ) and squares give E p . 
Dashed line is the Bogoliubov energy e p . Differences between 
diamonds and squares are due to non-Lorentzian structure in 
Q p . There are no free parameters in the theoretical results. 



For a fixed iVo we first solve the static GPE of Eq. © 
with n = m = / = 0to obtain < I ) o( r )- We then cal- 
culate and store the quasiparticle basis functions Uj(r) 
and Vi(r) and unperturbed energies q from the static 
limit of Eq. © for all states up to an energy cutoff 
Pcut ~ 130hu> r . Using these we can construct all static 
and dynamic terms, defined by sums and integrals over 
various functions of the quasiparticles. The numerical 
calculation is difficult because of the need to deal simul- 
taneously with the phonon character of low-energy states 
and significant single-particle effects. We therefore use 
an accurate Gaussian quadrature scheme together with a 
large value of E cut and a semi-classical approximation at 
high energy Q . The final results are converged to within 
5 x 10 _3 fiav. Further details are given in |l6j . 

Results for the m = 2 and m = modes are compared 
to experiment in Fig^ As can be seen, the theory pre- 
dicts a significant downwards shift for the m = 2 mode. 
The agreement with experiment is reasonable if we con- 
sider the temperature error in the experiment (of order 
5 — 10%) which is not shown. The downward curvature of 
the results is due to the scaling of the temperature axis 
from absolute to reduced temperature. For k^T ;§> Ao 
the shift is linear in T, as expected theoretically 

If we neglect thermal driving then similar behaviour 
is seen for the m = mode, as found in previous calcu- 
lations H, Qj 13 ■ Including this effect gives very differ- 
ent results, however, and the theory correctly reproduces 
the sharp upward shift in the excitation energy around 
t = 0.6. This is because an r 2 perturbation couples 
strongly to single-particle modes with frequency differ- 
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FIG. 2: Modulus squared of the dimensionless resolvent FIG. 3: Theoretical decay rates (T) compared with experi- 
H p (uj) x hu r as a function of frequency for t = (solid), men t for (a) m p = and (b) m p = 2. Symbols as in Fig.Q 
t = 0.66 (dot-dashed, x4) and t = 0.89 (dashed, x4). 



ences of 2uj r so the non-condensate response is peaked 
in this region. The effect on the condensate can be seen 
by plotting the modified resolvent 1Z as a function of fre- 
quency and temperature as in Fig. [21 The appearance 
of a growing peak at ui = 2uo r is due to direct driving of 
the non-condensate and is absent in an equivalent plot 
of Q. In this case the perturbation mainly excites the 
non-condensate which then drives the condensate, rather 
than the reverse. This explanation of the experimental 
results is consistent with the conclusions of BKS and JZ 
0, 0| . If the condensate drives the non-condensate the 
two oscillate out-of-phase. However, at high t the non- 
condensate drives the condensate above its resonance fre- 
quency and hence the oscillations are in-phase. The out- 
of-phase branch should be observable using a perturba- 
tion localized around the condensate. 

Fig. shows the results for the damping rates. Over- 
all the agreement with experiment is good, although the 
theory overestimates the damp ing rate at low tempera- 
tures. This was also seen by JZ |l(J and is possibly due to 
experimental difficulties in determining the temperature 
when the non-condensate fraction is small. For m = 0, 
the damping rate is underestimated at high temperature 
if direct driving of the thermal cloud is included for rea- 
sons which are currently unclear. 

In conclusion, we have presented numerical results 
from a gapless theory of condensate excitations which 
includes the anomalous average, Beliaev and Landau pro- 
cesses and the dynamic coupling of condensate and non- 
condensate fluctuations. Good agreement with the JILA 
experiment Q is found for the energies and decay rates 
of the lowest- lying states with m = 2 and m = 0. This 
shows that a consistent perturbative approach is capable 
of explaining the experimental results, contrary to state- 
ments in the literature LD, |l£J| ■ The anomalous behaviour 



of the m — mode is the result of direct excitation of 
the non-condensate by the external perturbation. 
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